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ABSTRACT 



CLc We present the results of three-dimensional simulations of the deep convective 

O ' envelope of a young (10 Myr) one-solar-mass star, obtained with the Anelastic 

■ Spherical Harmonic code. Since young stars are known to be faster rotators than 

, i^, ' their main sequence counterparts, we have systematically studied the impact of 

the stellar rotation speed, by considering stars spinning up to five times as fast 
^ ■ as the Sun. The aim of these nonlinear models is to understand the complex 

^ . interactions between convection and rotation. We discuss the influence of the 

Q\ • turbulence level and of the rotation rate on the intensity and the topology of 

the mean flows. For all of the computed models, we find a solar-type superficial 
differential rotation, with an equatorial acceleration, and meridional circulation 



. that exhibits a multicellular structure. Even if the differential rotation contrast 

Afl decreases only marginally for high rotation rates, the meridional circulation 



intensity clearly weakens according to our simulations. We have also shown that, 
for Taylor numbers above a certain threshold (T^ > 10^), the convection can 
develop a vacillating behavior. Since simulations with high turbulence levels 
and rotation rates exhibit strongly cylindrical internal rotation profiles, we have 
considered the influence of baroclinic effects at the base of the convective envelope 
of these young Suns, to see whether such effect can modify the otherwise near 
cylindrical profiles to produce more conical, solar-like profiles. 
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Introduction 



Young stars are the subjects of intense observational and theoretical research. There is 
a diversity of interesting objects, from massive Herbig stars to low-mass T Tauri stars. We 
focus in this paper on low-mass stars. During the pre-main sequence (PMs) they evolve from 
the birthline along the Hayashi track; during this phase they are fully convective. Next a 
radiative core develops and grows until the stars reach the zero-age main sequence (zAMs) 
and start burning Hydrogen. Thus energy transfer by convection takes place in a large 
portion of the PMS stars and, as such, is a key process for understanding their structure and 
their evolution. 

During their youth, these stars are surrounded by a disk, more or less dense and thick 
according to their evolution stage, with possible accretion from the disk to the star. Stars 
and disks are generally magnetically coupled. Because of the nature of these interactions, 
there is a large range of stellar rotation speed. There are slow rotators — spinning at the 
solar rate — and stars spinning a hundred times faster. As a direct consequence, stars on 
the ZAMS exhibit a very large dispersion of equatorial velocities, but T Tauri stars present 
more moderate rotation periods (typically between 2 and 15 days). The velocity distribution 
can be well modeled by considering the interaction between the star and its surrounding disk 
(IBouvier et al.l 119971 ). The rotation speed for a star on the ZAMS directly depends on the 
duration over which the central body and the circumstellar material have been coupled in 
earlier phases. During the main sequence, because of angular momentum losse s through a 
stella r wind, stars slow down and their rotation speeds follow Skumanich's law (jSkumanich 
19721 ). Thus young stars are generally faster rotators than their more evolved brethren. 
Such higher rotation rates should impact the internal dynamics of the stars and this impact 
must be studied with care. In order to progress on our understanding of such objects, one 
has to rely on both precise modeling of their internal structure and dynamics, and accurate 
observations. 



Three-dimensional hydrodynamic simulations are useful tools to understand complex 
and intricate physical processes, like the interaction between turbulent convection, fast 
rotation and possibly magnetic fields. Successes of such simulations in reproducing so- 
l ar observations, especially internal constraints on rotation provided by helioseismology 



fMiesch et al.ll2000l: Elliott et al.ll2000l : iBrun fc Toomrell2002l : iBrun et al.l 12004 : iMiesch et al. 



20061 ) . encourage us to pursue this effort in modeling other stars. 3-D simulations of convec- 
tive cores of A-type stars have thus been performed. The properties of convection and 



angular momentum redistri 



jution were analyzed in purely hydrodynamical models first 



(IBrowning. Brun. fc Toomre 



20041), then dynamo processes were studied with MHD simula- 



tions (IBrun. Browning, fc Toomrdl2005l ). Simulations of fully convective stars, like T Tauri 
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ones, have also been performed by iDobler. Stix. fc Brandenburg (120061 ). They have per- 
formed 3-D MHD simulations of a convective sphere embedded in a Cartesian box, aiming 
to model dynamo processes in fully convective stars. In this paper, our main objective is to 
continue this effort by performing 3-D hydrodynamical simulations of young solar-like stars. 

Today's observations are mostly limited to the stellar surface, even though with the re- 
cent launch of CoRoT, deep probing of the internal structure will be soon accessible. Despite 
being currently unable to reach the accuracy of observations of the solar interior given by he- 

lioseismic inversions, we can expect that asteroseismo logy will provide constraints on physical 

processes, such as the convection in the envelopes of stars (jMonteiro. Christensen-Dalseaard. &: Thompson 



2000| : iMazumdar fc Antial[200ll : [Roxburgh fc Vorontsovll200ll:lBallot. Turck-Chieze. fc Garcia 



20041 ). especially the young stars in open clusters (jPiau. Ballot. &: Turck-Chiezell2005l ). Among 
the current available proxies used to infer the inner part of young stars. Lithium is playing 
a central role. Lithium is burned by proton capture at temperature around 2-3x10^ K. 
Thus, its surface abundance is very sensitive to the temperature at the base of the convec- 
tive zone (bcz) and the physical p rocesses — diffusion, mixing — occurring in this region 
(iBrun. Turck-Chieze. fc Zahnlll999l ). It provides an efficient mean for checking deficiencies 
in classical modeling of PMS-stell ar structure. Indeed the evolu tion of ^Li abundance during 
the PMS is not so easy to predict. iPiau &: Turck-Chiezd (120021 ) have shown that ^Li burning 
depends strongly on many ingredients during its short phase of burning — typically from 2 to 
10 Myr for young suns — which corresponds to the contraction phase and a large increase of 
the temperature and density at the BCZ. Models of solar-like stars show a stronger d epletion 
of Lithium than the observed ones ( Ventura et al. IQQsi ^ Piau &: Turck-Chieze 2002), with a 
great dependence of the internal composition (see also Sestito et al. 20061 ) and of dynamical 
processes. The convective zone is clearly a crucial element and many questions have been 
addressed on the role of magnetic field on its extension and also of the effects of turbulence 
in this region. In this context, it is time to perform 3-D simulations of young stars, and see 
what we can learn from these models, first when convection and rotation interact, then when 
magnetic field is also considered (upcoming paper). 

In order to constrain these 3-D simulations, we summarize what is our current knowledge 
on their dynamics. The rate of differential rotation in latitude is known for many stars. Sev- 
eral techniques exist to infer this differential rotation; they mainly rely on the stellar magnetic 
activity. The first method is to use the migration during activity cycles of spots which rise 
through the surface at higher latitudes at the beginning of a cycle than at the end, according 
to the solar paradigm. Observations of luminosity fiuctuations indicate the rotation period 
for the mean latitude where spots emerge. Thus studying the period modulation within 
cycles gives a mea surement of differential rotation. Long term photometric monitoring by 
Henry et al.l (119951 ) revealed such secular modulation in period signature. They deduced that 
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the contrast of differentia l rotation AQ is almost independent of rotation {AQ oc 
Messina fc Guinanl (120031 ) have also deduced from long-term photometry the differential ro- 
tation rate for few young solar-like stars and found a correlation AQ oc J^o.ssio.s jg 
marginally significant relative to the error bars. Using Mount Wilson Survey of chromo- 
spheric Ca II emission-line fluxes for an hundred of stars, iDonahue. Saar. fc BaliunasI (119961 ) 
have found a strong dependence AQ oc Q^-'^^^-^. By using another metho d which consists of 



l ookin g for differential rotation as variations of stellar absorption lines, iReiners fc Schmitt 
(I2OO3I ) also have obtained a positive correlation AQ oc f20-53±o.35_ 



However, by focusing on young solar-type stars, ICameron et al.l (1200 ll ) have found a 
relation AQ ~ constant. This relation is obtained by processing separately K-type and G- 
type stars. AQ does not vary signifi cantly inside e ach c lass but is larger for G stars than for 
K ones. A more systematic study of iBarnes et al.l (120051 ) . indicates such a strong dependence 
on effective temperature but a marginal sensitivity to rotation {AQ oc Q^'^^^^'^^) . This study 
is mainly based on another, more direct, technique, wh ich consists of following displacem ents 



isp 

of spots on the stellar surface using Doppler imaging (iPetit. Donati. fc CameronI l2002h . 



Up to now, theoretical predictions of differential rotation in convectiye envelopes are 
maiii l y produced from mean field models (IKitchatinov &: Riidige ] |l995l : ll999l:lKuker Sz Riidigen 



1997 



Riidiger et al.lll998l : iKiiker fc Stixll200ll : iRempel 



2005a 



a) 



In this paper we aim to analyze the dynamics in the thick convective shell of a young 
(10 Myr) solar-type star, and particularly the influence of rotation on the convective motions 
and resulting mean flows — differential rotation, meridional circulation — in such a star, in 
order to provide scaling laws useful for observers and 1-D stellar evolutionary models. After 
defining the equations solved by the anelastic spherical harmonic (ASH) code and giving 
specifications of our different models (§[2]), we describe the spatial and temporal properties 
of convection in these simulations (§ Differential rotation profiles and their evolution 
with Q is discussed in § H] while § |5] is dedicated to the meridional circulation. The physical 
interpretation of all these results is discussed in § [61 Section [7] is focused on a peculiar form 
of vacillating convection we have found during this study. In conclusion (§ [8]) main results 
are summarized and perspectives are proposed. 



2. Posing the Problem 

2.1. Anelastic equations 

The simulations described here were performed using the Anelastic Spherical Harmonic 
(ASH) code. ASH solves the three-dimensional anelastic equations of motion in a rotat- 
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Miesch et al. 


2000; 


Brun et al. 


2004) 



ables and linearized in thermodynamic variables with respect to a spherically symmetric 
mean state. This mean state is taken to have density p, pressure P, temperature T, specific 
entropy S; perturbations about this mean state are written as p, P, T, and S. Conservation 
of mass, momentum, and energy in this rotating reference frame are therefore written as 

V ■ ipv) = 0, (1) 

p{^ + (vV)v + 2n,xv)=-VP 

+pg-V-V-[VP-pg], ^ ' 

pf§ = V ■ l^rPCpVif + T) + KpfV{S + S)] 

-pfv-V{S + S) (3) 
+2pu [eijCij - 1/3(V ■ vf] , 

where cp is the specific heat at constant pressure, v = {vr,ve,v^) is the local velocity in 
spherical geometry in the rotating frame of constant angular velocity fto, 9 is the gravita- 
tional acceleration, k,. is the radiative diffusivity, and T> is the viscous stress tensor, with 
components 

Vij = -2pu[ei^ - 1/3(V ■ v)6i,], (4) 

where 6,^ is the strain rate tensor. Here u and k are effective eddy diffusivities for vortic- 
ity and entropy. To close the set of equations, linearized relations for the thermodynamic 
fiuctuations are taken as 

p P T P S 



assuming the ideal gas law 



(5) 

p P T -fP Cp' ^ ^ 



p = npf, (6) 



where TZ is the gas constant. The effects of compressibility on the convection are taken into 
account by means of the anelastic approximation, which filters out sound waves that would 
otherwise severely limit the time steps allowed by the simulation. 

Convection in stellar environments occurs over a large range of scales. Numerical simu- 
lations cannot, with present computing technology, consider all these scales simultaneously. 
We therefore seek to resolve the largest scales of the nonlinear fiow, which we think are 
likely to be the dominant players in establishing differential rotation and other mean prop- 
erties of the convection zone. We do so within a large-eddy simulation formulation, which 
explicitly follows larger scale flows while employing subgrid-scale descriptions for the effects 
of the unresolved motions. Here, those unresolved motions are treated as enhancements to 
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the viscosity and thermal diffusivity (z/ and k), which are thus effective eddy viscosities and 
diffusivities. For simphcity, we have taken these to be functions of radius alone, and to scale 
as the inverse of the square root of the mean density. We emphasize that currently tractable 
simulations are still many orders of magnitude away in parameter space from the highly 
turbulent conditions likely to be found in real stellar convection zones. These large-eddy 
simulations should therefore be viewed only as indicators of the properties of the real flows. 
We are encouraged, however, by the success that similar simulations (cf. §[1]) have enjoyed in 
matching the detailed observational constraints for the differential rotation within the solar 
convection zone provided by helioseismology. 



2.2. Numerical approach 

Velocity and thermodynamic variables within ASH are expanded in spherical harmonics 
Y^{6,(j)) in the horizontal directions and in Chebyshev polynomials T„(r) in the radial. 
Spatial resolution is thus uniform everywhere on a sphere when a complete set of spherical 
harmonics of degree i is used, retaining all azimuthal orders m. We truncate our expansion 
at degree £ = £max, which is related to the number of latitudinal mesh points Ng as £max = 
{2Ng — l)/3, take A'"^ = 2N0 latitudinal mesh points, and utilize Nr collocation points 
for the projection onto the Chebyshev polynomials. The grid resolution Nr x Ng x N^ 
(cf. Table [I]) we have considered depends on the degree of turbulence of the model. An 
implicit, second-order Crank-Nicholson scheme is used in determining the time evolution of 
the linear terms, whereas an explicit second-order Adams-Bashforth scheme is employed for 
the advective and Coriolis terms. The ASH code has been optimized to run efficiently on 
parallel supercomputers such as the HP ES45, on which the simulations of this work have 
been performed. 



2.3. Setting the young sun model 

We describe and analyze in this paper a simplified 3D model of the convective envelope 
of a young one-solar-mass star where the plasma is composed only of a mixture of fully 
ionized Hydrogen and Helium. This model is built from the structure of an accurate 1-D 
stellar model (model Y) to get realistic values for the radiative opacity, density, and pressure 
profiles. 



Model Y is obtained with the CESAM stellar evolution code (IMorell 119971 ) for an age 
of 10 Myr from the birthline. At this age, the star is located at the turning point of the 
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PMS evolutionary track, just after the Hayashi descent. Its radiative cor e is developing 
quickly and there is a strong lithium burning (jPiau &: Turck-Chiezd |2002| . Fig. 1). It is 



characterized by an effective temperature of 4590 K, a luminosity = 1.86 x 10 ergs 



(0.48 Lq) and a radius R^, = 7.67 x lO^'' cm (1.1 R©). The outer convective zone covers 
45% of the stellar radius or 37% of the total mass. Such a convective envelope is almost 
twice as thick as the solar one, and, more importantly, 15 times more mas sive. For such a 
model, we use t he OPAL equation of sta te (IRogers et al.lll996l : lRogersll2000l ) and the nuclear 



Jl998h ;o describe the microscopic properties of the stellar 



cross sections of lAdelberger et al. 

plasma. OPAL opacities (llglesias &: RogersI Il996l ) are completed at low temperature by 



tables of [Alexander fc Ferguson! (119941 ). Convection is computed using the classical mixing 
length treatment calibrated to match the observations of the Sun at 4.6 Gyr [a ~ 1.9). 

The simulations performed with the ASH code were initialized from this model Y. The 
gravity g, radiative diffusivity Kr, and mean density p radial profiles are the starting points 
for an iterative Newton-Raphson solution of the hydrostatic balance and for determining 
the gradients of the thermodynamic variables. The mean temperature T is then deduced 
from equation [6l This technique yields mean profiles in good agreement with the initial 1-D 
stellar model Y (Fig. [T]). 

The computational domain extends from about 0.55 to 0.95 i?*. The overall density 
contrast in radius is around 60, implying to noticeable compressibility effects. The simula- 
tions do not model the outer most layers of the star, where the convective scales become very 
small due to the sudden decrease of the density scale height. This choice excludes the H and 
He recombination zone and the superadiabatic layer for which the microscopic description 
is more complex. The inner boundary corresponds to the base of the convective zone; the 
presence of a potential overshooting layer is not taken into account in this work. 

We impose impenetrable and stress-free conditions for the velocity field and a constant 
flux — i.e., a constant entropy gradient — at both inner and outer boundaries, except for 
the model Yc5S discussed in §[71 In this model, an entropy gradien t in latitude is imposed at 
the bo ttom of the shell to model the influence of a tachocline as in iMiesch. Brun. &: Toomre 



2.4. Early phases and relaxation of 3-D models 

Some simulations have been started from a quiescent state with a uniform rotation; 
others use evolved solutions in which we modify diffusivities or rotation rate. One can track 
the development of the convection by following the mean enthalpy flux crossing the shell. 
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The enthalpy flux is defined as 

Fen = VrpCpT (7) 

Figure [2] shows the temporal evolution of Fen, averaged over the full domain, along with the 
evolution of mean kinetic energy KE, and its axisymmetric and non-axisymmetric compo- 
nents (see § 16.11 for more details). Given the high Rayleigh numbers used in all our models, 
convective instabilities develop from the quiescent state and grow exponentially over about 
500 days, until saturation. After a phase of ~1500 days over which relaxation oscillations 
occur, statistically stationary states are usually reached. Only statistical fiuctuations modu- 
late the temporal trace in the late evolution. When such a state has been attained, we have 
integrated the model about a thousand days, namely several convective-overturning times 
and several rotation periods, even for the most consuming simulations, for a total of one-half 
million node hours spent. 

One can see for example how the energy is transfered from the bottom to the top of 
the shell, on average, for a relaxed simulation. Sev eral processes tak e part in this transfer 



(for detailed expressions, see equations [11]-[16] in iBrun et al.ll2004l ). Figure [3] shows the 
horizontally averaged distribution of each fiux in a typical model (Ycl), after averaging over 
time (1200 days). The main component is of course the enthalpy luminosity Le„ advected 
by convection (luminosities L are defined from fiuxes F by the relation L = Anr'^F). Close 
to the bottom, radiative luminosity Lrd becomes significant because of the steady increase 
of radiative conductivity Hr with depth. Near the top of the layer, the heat transport is 
dominated by the subgrid-scale turbulence that yields the term Led- This fiux is proportional 
to a specified eddy diffusivity and to the mean radial gradient of entropy and serves to carry 
the total fiux through the upper boundary and prevents the entropy gradient there from 
becoming superadiabatic compared to the scales of convection that we are able to resolve 
spatially. There are two minor extra contributions. The first is the kinetic energy fiux {F^e), 
which is slightly negative, except when Led dominates. This negative value i s a consequence 



of the fast downfiow sheets and plumes achieved by effects of compressibility (IHurlburt et al. 



19861 ). The second is the viscous flux Lj,, which contributes to the total flux mainly near the 



top where there is a viscous layer. Proflles are very similar in the other simulations. 

In this work we will concentrate on three series of models, respectively called Ya, Yb 
and Yc. These models possess different levels of turbulence and Prandtl numbers Pr=4, 1 
or 1/4 {Pr = v/k). For each Prandtl number considered, we have performed models with 
f2o = 1, 2 and 5f20 {VLq = 2.6 x 10~^ rads~^, corresponds to a period of 28 days), for a 
total of nine models. We have also computed flve supplementary models spinning at 5 ^0 , 
with higher turbulence levels, different Prandtl numbers or boundary conditions. This study 
spans a large range of Taylor and Rayleigh numbers, Ta G [10^, 10^°] and Ra G [10^, 10^]. All 
the model parameters and characteristics are detailed in Table [H 
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3. Convective patterns and their evolution in deep spherical shell 

Convection in young solar stars is a major player since it transports heat outward from 
the deep internal parts, redistributes angular momentum and establishes large scale mean 
flows. Understanding these intricate physical processes is thus crucial if we want to progress 
in our knowledge of such stars. We here start by describing in details the convective patterns 
established in our simulations, how they are modifled by the rotation rate and the turbulence 
level and how these nonlinear interactions lead to complex temporal evolutions. 



3.1. Topology of the convection 

Figure m shows maps of radial velocity v^. and temperature fluctuation T near the surface 
and at mid-depth of the convective envelope for models of the Ya series. Near the surface, 
the velocity fleld of the Yal simulation exhibits networks of downward flows conflned to the 
periphery of convective cells, and broader and weaker upward flows in their centers. This 
asymmetry is made possible by the anelastic approach which includes compressibility effects. 
This quite laminar model has a "classic" structure composed, in the equatorial area, of so- 
called ban^j]^^_celkj_aligned with rotation axis. It is a well-known pattern in solar modeling 
(e.g., 



Miesch et al.l 120001 ). Around the poles, convective cells look less elongated. 



By comparing maps at different depths, one recovers some patterns at similar places, 
which shows the great connectivity of the downflow network across the whole convective 
domain. In fact, upward and downward flows are not radial, i.e. aligned with the local 
gravity, but are tilted by the rotation and tend to be aligned with the polar axis. In Yal 
model, upward structures can start, in a meridional plane, at the equator at the bottom of 
the shell and reach the top near a latitude of 30°. Around the poles, tilts of structures are 
less signiflcant because Q,o are g are almost aligned. 

Comparing radial velocity Vr and temperature fluctuation T maps shows their strong 
correlation. It just means that hotter material goes up and cooler material sinks. This strong 
correlation of Vr and T leads to large outward enthalpy flux and implies heat transport by 
convection (see Fig. [3]). 

Temperature maps show also a clear variation with latitude. This zonal structure varies 
with depth. At mid-depth (as at the bottom, not shown), there is a monotonic variation from 
the hotter poles to the cooler equator. Close to the surface, the temperature increases again 
in the equatorial region. The amplitude of these variations is around 2K. This axisymmetric 
thermal variation between poles and equator plays an important role in the understanding 
of differential rotation proflles (see discussion §[6]). 
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When the rotation speed is increased (see case Ya5), the pattern shape is modified 
but retains its overall organization. Asymmetry of fiow is reduced because the convection 
is less vigorous. Indeed, we notice that speeds of convec tive fiows have been significantly 



reduced, because rotation tends to stabilize fiows (e.g., IChandrasekharl Il96ll ). Stronger 
spatial variations in fiow velocities also appear, which can be seen on the Vr maps as variations 
on the image contrast. Another noticeable difference is the pattern size. In the Ya5 model, 
convective patterns are clearly smaller than in Yal. This is also a well-known effect, which 
can be understood simply by following the linear growth of co nvective instabilities: degree s 



£c of the most instable modes increase with the Taylor number (iGlatzmaier fc Gilmanll 19811 ). 
From the maps, we can estimate that ^c ~ 18 in Yal and increases to 42 in Ya5. We have 
performed an extra model thus Ya5T, similar to Ya5 but more turbulent of which Ta is 16 
time higher; in this case, estimated £c goes up to 60, producing even finer structures. In 
addition, an important impact of rotation on fiows is the tilts of structures are strongly 
pronounced in the equatorial area. Convective structures are essentially aligned with the 
rotation axis in the Ya5 case (and somewhat less in Ya5T). In this case, the Coriolis term 
dominates the o thers; convec tion tends to develop along Busse's columns as a consequence 



(see for example iBussd |2002| ) . 



Finally, the rotation does not affect the structure of the temperature fluctuations. The 
zonal structure of T is very similar to those of Yal in that, the same kind of bands are 
visible; only the thermal contrasts have been reinforced, especially for Ya5T, in agreement 
with the stronger zonal flow present in this case (see §|l]and[n]). 

Convective motions obviously modify the radial structure. In a deep convective zone, 
the radial entropy gradient is negative (unstable region), but close to zero, especially at the 
bottom of the shell where the convection layer is close to adiabatic. Table [1] shows that 
for models with exactly the same diffusivities v and k, the Rayleigh number increases with 
rotation rate. It indicates that rotation modifles dS/dr; that is confirmed by Fig. El which 
shows the entropy gradient for the Ya model series. As described above for the velocity fields, 
the stabilizing effects of rotation reduce the vigor of the flow. As convection is less vigorous, 
it becomes less efficient and the entropy gradient needs to become steeper to evacuate energy. 
For the same reason, since the intensity of convection is greater in the simulation Ya5T, more 
turbulent than Ya5, the entropy gradient is less strong in Ya5T. 

We now compare the Ya series with the series of simulations Yc to see the effects of 
decreasing the Prandtl number from = 4 to = 1/4. Convective patterns of the three 
main Yc simulations are shown Fig. O analog of Fig. |H In going from the Yal to the Ycl 
case, the Taylor number is increased, as it was when Qo was increased. However, in this case 
there is no visible stabilization effect since, thanks to the reduction of z/, the turbulence level 
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increases as well. Ycl is the most turbulent run that we performed, if we exclude the special 
cases Yc5T and Yd5, as shown by the Reynolds numbers Re (Tabled]). Banana cells are less 
visible and structures are more complicated. Relative to Yal, the asymmetry of up/down 
flows is more pronounced in Ycl. The small-scale vorticity is more intense in this model than 
in more laminar ones, as suggested by a higher Rossby number Rq (Tabled]). By comparing, 
on Vr maps, the convective patterns in the polar and equatorial regions, it becomes clear 
in this simulation that the behaviors of convection in these regions are different. We can 
roughly separate what occurs inside and outside the cylinder tangent to the inner shell and 
parallel to the rotation axis. Inside this cylinder, convection follows a "polar regime" slightly 
influenced by rotation, outside it follows an "equatorial regime" where rotation influences 
the convection much more. 

When the rotation rate VLo is increased, convective features are modified. We can see that 
in the equatorial region convection tends to be localized, favoring a limited band of longitudes 
at a given time. By contrast, at the solar rotation rate, this is not very pronounced. We note 
that around the equator there is a zone (close to the map center) where the velocity is slightly 
stronger. When the rotation is double or more (Yc2, Yc5), vigorous convective motions are 
concentrated in a small region in longitude. In the rest of the equatorial zone, convective 
velocities are markedly smaller (by a factor 5 or 10 typically). Polar convection is much less 
affected. We observe that, in this region, convective patterns are smaller, for reasons we 
have previously developed. Thus, convective motions inside and outside the vertical tangent 
cylinder are only loosely connected. 

Finally, we can compare once again Vr and T maps in Fig. [51 Correlations are apparently 
less obvious. This is because the axisymmetric component of T presents a stronger contrast 
than for Yal or Ya5 (it is quite similar to Ya5T). If the m = components are removed from 
T, correlations clearly appear again. We retrieve in these model runs the same structure of 
azimuthal bands as in the previous series. This is a characteristic common to all our model 
simulations. Only the thermal contrast of these bands varies. It as been increased here by 
a factor 4~5 relative to Ya series. As for the previous series, the thermal contrast increases 
with VLf,. 



3.2. Pattern temporal evolution 

After having described the spatial structures of convection in our simulations, we discuss 
here how they evolve with time. To do so, we have followed the evolution of convective 
patterns thanks to a time-longitude (0,t) diagram. Figure [7] shows such diagrams for Yal 
and Ycl models. For given radius and zenith angle 6'c, we have plotted the radial velocity 
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field Vr in the (<^, t) plane, serving to follow, for each longitude, the temporal evolution of Vr- 
We have also plotted "shifted" diagrams, in {(f)s,t) coordinates, with 



'c J to 

where v^i, is the azimuthal average of v^, in order to compensate for the local advection of 
the differential rotation. 

A first look at the non-shifted maps (Figs [7h and Uh) shows that velocity structures 
propagate eastwards. It is still moderate for Yal, but becomes really strong for Ycl, so 
much so that the map becomes difficult to read. For such cases it is very useful to consider 
shifted maps since they are more readable. Since the maps are plotted in the corotating 
frame, such an eastward propagation is a signature of a strong differential rotation, where 
the equator is accelerated relative to the mean rotation speed (cf. §|1]). 

Let us focus on Yal. Maps show convective patterns are persistent during several 
turnover times: we can follow them in the map during the complete time sequence. Figure [7^ 
shows that differential rotation does not account for the total eastward propagation. There 
is still a small residual drift. We thus see a self-propagation of convection structures. For 
example, the structure located at = —80° at the initial time is around —60° longitude 300 
days later. We clearly see also the regular appearance of downfiow lanes and their ongoing 
mergers. These maps show undoubtedly the temporal modulation of speed and thickness of 
flows, as a consequence of the intrinsic fluctuation of convective motions. 

It is instructive to make a direct comparison with the pattern evolution seen in Ycl 
simulation. The shifted map of this turbulent model (Fig. UU) reveals conspicuous temporal 
modulations, stronger than in Yal. Moreover, the persistence of downfiow lanes is shorter. 
The hot spots, described in previous section, are persistent patterns. On the shifted map, 
they obviously drift westward, a signature of retrograde motion. 

The localized convective structure seen on Yc2 velocity maps (Fig. [5]) is also persistent 
and presents the same kind of retrograde motion (not shown). In the Yc5 simulations, 
the localized convective patterns also similarly propagate, but their temporal evolution is 
very peculiar: the convective level fluctuates highly in this case, so much so the pattern 
periodically fully disappears, reappearing later. Section [7] is dedicated to the study of this 
special case. 
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4. Differential rotation produced 

We have seen in the previous section that convective patterns propagate in longitude 
under the influence of a large-scale mean longitudinal flow. Here, we discuss in more detail 
this mean flow: the differential rotation, and how it varies, both in contrast and proflle, with 
the rotation rate. 



4.1. Differential rotation profiles 

Figure [8] shows the temporal and longitudinal averages of angular velocities achieved for 
six representative models. To make comparison with observations easier, we have converted 
mean azimuthal velocities i)^ into sidereal angular velocities VL = VL^ + v^/ {r smO). The 
rotation proflle of our reference case, Yal (Fig. [Hh), exhibits a clear equatorial acceleration, 
similar to the Sun. All of our models present such a behavior. I t is a consequence of the 



domination of the Coriolis force over the buoyancy (IGilmanl 119771 ). as measured by the low 
convective Rossby number : Table [1] indicates that Roc < 1 in every case. In Yal, the 
contrast of this differential rotation is not very large, but the decrease of Vt from the equator 
to the poles is monotonic (see right panel). The rotation proflle is rather conical, i.e., radial 
cuts at different latitudes are independent of r. 

On Fig. [HI we can see the effects of increasing T^, following two different paths : flrst 
by decreasing the viscosity v and so Pr : simulations Yal to Ycl via Ybl (panels a, 6, c); 
second by increasing VLo in going from Ycl to Yc5 via Yc2 (panels c, rf, e). 

In going from Yal to Ybl, we flrst notice the increase of the contrast of the differential 
rotation (Fig. W?)-, multiplied roughly by a factor of 4. The proflle remains monotonic with 
respect to the latitude coordinate. However, it becomes less conical, and tends to becomes 
more cylindrical, i.e. with isorotation contours parallel to the stellar axis. Compared to Ybl, 
in Ycl (Fig. [HI:), the proflle becomes even more cylindrical and the contrast is increased by a 
factor of 2. Such a contrast of 130 nHz is typical of those observed at stellar surfaces (around 
90nHz for the Sun). The rotation rate decreases monotonically from equator up to 60°, but 
the monotonicity is lost in the polar caps. Moreover the behavior of the rotation is different 
in the northern and southern polar caps. As we have already seen in § [3] concerning the 
convection, dynamics inside and outside the tangent cylinder are practically disconnected, 
and the two polar regions are disconnected too. However, very long averages should result 
to identical proflles. 

By following the path #2 (Ycl, Yc2, Yc5), the cylindrical aspect is more and more 
marked and the isorotation contours for Yc5 model (Fig. [Ht) are almost purely vertical. 
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The rotation behavior in the polar regions varies from one model simulation to the next 
in sequence, from one hemisphere to the other. Nevertheless, outside the tangent cylinder, 
the equatorward monotonic increase of Q is always present. The contrast of the differential 
rotation is slightly reduced along this path. Even though it is not very pronounced for this 
Yc series, it is more noticeable for Ya and Yb model series. Thus, models with high Pr and 
high O reach a rather weak differential rotation. The extreme case, Ya5, rotates practically 
as a solid body. 

Along both paths, when Ta is increased, rotation profiles become more cylindrical. How- 
ever, the contrast of the differential rotation responds differently in the two cases: it rises 
along path ^1 and declines when increases. This is an effect of the stabilizing role of ro- 
tation, an effect already mentioned in the previous section. This means that a non negligible 
part of the properties of linear unstable modes is retained in these non-linear simulations. 

We will discussed of the rotation profile of Yc5S case and the temporal modulation of 
Yc5 profile further in § [71 

4.2. Rotation scaling law 

We have seen that differential rotation contrasts can vary with variations in the model 
parameters. To quantify these variations we have computed for each model the contrast AQ 
at the surface of the convective shell, between the equator and the latitude 60°(Figs. [9]& [T0|) . 
We have chosen to compute it near the surface because That is where rotation is observed 
in stars (cf. §[I]). Moreover, by choosing to cut the range at 60° we avoid the polar regions 
where peculiar behavior occur in Yc series. These contrasts AQ are summarized in Table [2] 
among other useful quantities. We have plotted AQ as a function of the turbulence level 
characterized by Reynolds number R'^ (Fig. [9]), and of the rotation rate Qo (Fig. [lO]) for all 
our models. Thus we get a global vision of the sensitivity of AQ to different parameters. 

First, there is a clear and strong dependence of AQ on the turbulence level. The 
simulation Ya5T is more turbulent than Ya5 — with the same Pr and Qq — and reaches a 
stronger contrast, even higher than Yal (see Table [2]). That is also true when we compare 
Yb5T to Yb5 or Yc5T to Yc5. To show the strong impact of turbulence level, we can 
compare AQ and R'^. We have defined two Reynolds numbers: R^ = vL/u and R'^ = v'L/u 
(see Tabled]). The former is computed with the rms velocity at mid-depth layer w; the latter 
is computed with v' , which is the rms velocity of convective motions, i.e. that evaluated after 
removing the azimuthal averages to the velocity components (Table [2]). It is more relevant to 
compare AVt with R'^ than with Rf.: when the differential rotation is strong, v^p is the major 
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contribution to v and so increases the value of i?e, thus, by construction, we have to find 
a correlation between and i?e- On the contrary, R'^ takes into account the convective 
motions only, since the differential rotation has been filtered for its evaluation. Thus by 
comparing Afi to -Rg, we really measure the impact of the convective turbulence level on 
the differential rotation. Figure [9] shows the clear direct correlation between Af2 and R'^ we 
have found. The contrast Af2 increases with R'^ and seems to saturate around 140 nHz in 
the more turbulent cases. The signification of this correlation will be discussed in § 16. 2[ 

The variation of Af2 by going from Ya to Yc, seen in Table [2], seems to indicate a 
dependence to Pr- However, this variation is mainly the effect of the growth of the turbulence 
level from Ya to Yc. Thus we have to compare models with same turbulence levels, such as 
Ya5T and Yb5T. The latter, with a lower P^, achieves a slightly higher differential rotation 
than the former. Nevertheless, by comparing Yc5T to Yd5, one could conclude the opposite. 
In any case, the sensitivity to Pr is obviously weak in comparison with R'^. 

Finally, we have carefully studied the effects of Vto on Af2. For each series Ya, Yb and Yc, 
we have deduced a scaling law from every set of three simulations having the same parameters 
except Vto (for instance the set Yal, Ya2, & Ya5). We have fitted a power law Af2 oc The 
resulting fits are plotted in Fig. [TOl The a exponent increases from around —0.5 for the more 
laminar simulations (Ya and Yb) to —0.19 for the more turbulent ones (Yc). Thus for every 
series Af2 decreases with the rotation rate, more or less strongly. When VLo is increased, due 
to stabilizing effects R'^ decreases and therefore AQ too. This phenomenon could explain the 
sign of the slope. However, the fall of R'^ is not the only effect: by comparing Ybl and Yb5T, 
we notice that, at 5f20, the turbulence level must be higher (i?g ~ 55 for Yb5T against 38 
for Ybl) to achieve the same contrast (Af2 k, 64 nHz). If we focus now on the Yc series, 
the decrease in this case is less pronounced. All of the models achieve contrasts between 100 
and 140 nHz, which is close to observations. Turbulence seems to be sufficiently developed to 
make the contrasts of this set of models less sensitive to R'^. The exponent of the scaling law 
we have obtained with this more reliable sequence is only slightly negat ive (—0.19). Despite 



differences in details, this is not far from the most recent results of iBarnes et al.l (120051 ) 
cited in the introduction. They have found a marginally positive coefficient which is mainly 
compatible with the relation Af2 ~ constant. If we extrapolate our sequence Ya, Yb, Yc to 
more turbulent series with lower P^, we expect that Af2 tends to be independent of Vto- 



5. Meridional circulation 



The meridional circulation is markedly smaller than the differential rotation, but can 
play a non-negligible role in the dynamics, especially in the heat and angular momentum 
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transport (see § I6.2p . Further, its topology is an import ant parameter for na odels of stellar 



dynamos, such as the Babcock-Leighton type ones (e.g.. |Jouve &: Brunll2007l ). 



Figure [TT] shows the time-averaged meridional circulations produced in three of our 
model simulations: Yal, Ycl and Yc5. We have plotte d in the meridional plane the mass 



flux streamfunction, as defined in iMiesch et al.l (j2000|, eq. [7]). We have also plotted on 



this figure the velocity of the circulation at the top of the domain. Since our boundary is 
impenetrable, the radial component Vr vanishes, so vq is the total speed of the meridional 
circulation. The amplitude of the flow is a fewms"^, ranging between 1 and 10 for the set 
of simulations considered in this study. Due to the compressibility effects, deep inside the 
convecting shell, velocities decrease strongly to conserve the mass flux in denser regions. 

Considering first the reference case Yal, we see a well-organized multicellular circulation, 
where both hemispheres present almost antisymmetric cells. If we exclude what happens 
around the poles, there are four cells along the radius. The cells near the top of domain induce 
poleward flows at the surface near the equator, that are similar to the solar observations and 
simulations. Near the poles, cells are smaller, vertically oriented, and go across the entire 
thickness of the modeled shell. 

Turning now to Ycl model with a higher Tq, we notice a more intricate and complex 
topology with an increase of the number of cells. The cells near the surface in the equatorial 
region retain the same circulation direction and a similar intensity, even if they are a little 
less extended. We notice that meridional circulation becomes stronger along the tangent 
cylinder. In this region, the structures become more "concentrated" and long and thin cells 
appear. If we continue to increase by increasing Qq, this effect is reinforced; long cells near 
the tangent cylinder are more concentrated and dominate in intensity over others further 
away from the tangent cylinder. We can find some similarities between this phenomenon 



and the Stewartson layers seen in simulations of rotating stellar radiative regions (IRieutord 



20061 ) . In the equatorial region, the cells already seen both for Yal and Ycl, are still present, 



having greater extent in latitude, but weaker velocity. 

To show the dependence of the meridional circulation on the rotation rate, we have 
scale it by the maximum velocity of the equatorial cells, present in all models. We plot the 
maximum of ve at the surface in these cells as a function of flo in Fig- [ISl This velocity 
is generally the highest vq produced at the surface of our models, except for Ycl (but that 
does not change our conclusion). We see that ■Ogmax in the Yb series are higher than in both 
Ya and Yc for rotations 1 and 2 times the solar rotation rate. This hierarchy between series 
is different from the one observed for AQ. However, if we should have considered the total 
kinetic energy in the meridional circulation (see MCKE in Table [2]), the hierarchy between 
Yb and Yc should be the same as for the differential rotation. The difference has been 
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mentioned before: in Yc series, a large part of MCKE is concentrated deep inside, along the 
tangent cylinder, and is not visible in ■Ogmax- 

As previously done with AQ, power laws are fitted to the three sets of model simulations 
with identical input parameters, excepting Qo- For all series, the power-law exponent is 
clearly negative (between —0.5 and —1). In contrast to what we have seen for differential 
rotation, this trend of declining "Oemax with Qo is rather strong. Even if in Yc5T, the more 
turbulent counterpart of Yc5, the intensity of meridional circulation slightly increases relative 
to Yc5, but it stays at a lower level than both Ycl and Yc2, a property that was not found 
for AQ. Thus, if the differential rotation decreases only marginally for high rotation speeds, 
the meridional circulation clearly weakens according to our simulations. 



6. Interpreting the dynamics 

Our simplified model of a rapidly rotating convective envelope yields a nonlinear dynam- 
ical system with complex feedbacks. Numerical simulations are a very efficient tool to assess 
the dynamical balances that are achieved in such systems. They provide detailed insights 
in the way kinetic energy, heat and angular momentum are exchanged among the different 
reservoirs. Given the importance of understanding the internal differential rotation profile of 
stars, it is necessary to identify the processes involved in establishing these large scale fiows. 
In the following we discuss how the differential rotation and meridional circulation discussed 
in § m and § [5] are produced. 



6.1. Energy balance 

We start by analyzing the energy balance of the various simulations discussed in the 
paper. In Table [2] we have summarized key quantities such as rms velocity or volume av- 
eraged kinetic energy (KE = |p(f^ + + v^)) achieved in the models. In order to have a 
more precise analysis, we decompose the rms velocities into their three components, with 
or without subtracting the azi muthal averages, a nd K E into its axisymmetric and non- 



axisymmetric components, as in iBrun &: Toomrd (l2002l ). We define the kinetic energy in 



differential rotation as DRKE = \pv'^, the kinetic energy in the meridional circulation as 
MCKE = ip({;2 + vj), and the kinetic energy in non-axisymmetric (convective) motions as 
CKE = KE - DRKE - MCKE. 

Turning first to the rms velocities we notice a large range of values, from a few up 
to about one hundred ms~^, especially for v and f;,^. Both these velocities have the same 



- 18 - 



magnitude in every model; in contrast, when comparing the longitudinal velocity v^p to its 
radial and latitudinal counterparts, we see significant differences: can be up to 10 times 
higher than Vr and vq. This is mostly due to the large scale differential rotation achieved 
in almost all the simulations and most certainly in the Yc and Y#T series as seen in § |H 
If we instead compare the rms velocities with their axisymmetric mean subtracted, we find 
that the three components have about the same amplitude with v[. ~ Vr being the fastest. 
This yields convective fiows that are largely isotropic. If we define the isotropy index as 
Qr = v[.'^/v''^, we find that it varies from about 1/2 in the Ya series to about 1/3 in the Yc 
series. At large Prandtl number the anisotropy of the fiow is thus greater than when Pr is 
small. This can perhaps be explained by the fact that the fiows are more turbulent in the 
latter series, having both larger Reynolds and Rayleigh numbers (see Tabled]). This leads to 
convection that is less dominated by large scale rolls (consistent with the disappearance of 
banana cell structures) that are characterized by large coherent radial motions. Instead, the 
convection, characterized by more chaotic random motions, is dominated by strong downward 
convective plumes as discussed in §[3l Within one of these plumes downwards vortical motion 
dominates, but when forming horizontal averages their small filling factor do not contribute 
as much as large scale coherent rolls. By increasing Qo, CLr tends to rise slightly, probably for 
similar reasons. The stabilizing effects of rotation also reduce the supercriticality of fiows. 
In this case, convection is less vigorous and more infiuenced by large scale rolls. 

Studying the partition of kinetic energy between its axisymmetric and non-axisymmetric 
components is also very instructive. As was shown in Fig. [2], our simulations, after a short 
linear convective instability phase, undergo a nonlinear saturation and finally reach, in most 
of the cases after about 2000 days, a statistically stationary state, over which meaningful 
temporal and volume averages of the kinetic energy and its mean and fiuctuating contribu- 
tions can be evaluated. Turning again to Table [21 we notice that the energy in the meridional 
fiows is very weak (less than 0.2% of KE) compared to DRKE and CKE that contain most 
of the energy. Again, this is because the meridional fiows in such systems result from small 
i mbalances be tween several large terms and are less easily developed than azimuthal motions 



(lMieschll2005l ). Another obvious trend, present in all the series Ya, Yb and Yc is that an 
increase of the rotation rate leads to a smaller kinetic energy available for the system. This 
reduction is about a factor of 3 for Ya and Yb series and only 1.6 for the Yc series. Rotation 
thus stabilizes convection and leads to weaker mean or fiuctuating fiows, as the decrease in 
absolute value of DRKE, CKE and MCKE indicates (note however that in terms of percent- 
ages CKE & DRKE remains about the same). The fact that the simulation in the Yc series 
are less sensitive to the increase of Qq is certainly due to their higher degree of supercritical- 
ity. We also find that the ratio of fiuctuating and mean (mostly azimuthal) motions changes 
in favors of the latter as we decrease the Prandtl number. In the Ya series with Pr = 4, 
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DRKE represents only about 35% of KE, with most of the energy concentrated in CKE. By 
contrast DRKE represents respectively about 80% and 95% for Yb and Yc series. It is thus 
clear that the higher is the Taylor number, the stronger is DRKE and the associated differ- 
ential rotation, as the column Af2top also clearly indicates (see also discussion in §|1]). Given 
the strong stabilizing effect of Qo on the motions, we have also performed simulations with 
rotation set at 5 times the solar rate, but with a higher degree of surpercriticality (models 
labeled with T). All these models posses a large DRKE, including model Ya5T. Thus the 
fact that DRKE is found to increase with decreasing is perhaps more linked to the super- 
criticality of the models than to the ratio between viscous and thermal diffusivities. A useful 
number to assess the influence of rotation on convection is the convective Rossby number 
Roc (see Tabled]). This number does not include u or n. We find that this number is always 
relatively small in the simulations with the largest differential rotation contrast, indicating 
a dominant effect of rotation over convection, meaning that as the convective flows reverse 
direction they are significantly tilted by the Coriolis force. 

In the Yc series, for the models rotating at 5 times the solar rate, the convection develops 
a rather intermittent behavior, in both space and time. Computing time-averaged quantities 
such as KE and rms velocities can thus be misleading, since they depend strongly on the 
state of the convection (quiescent or excited). In cases like Yc5 and Yc5T, we find that at 
the peak of the convection burst CKE can account for almost 30% of KE before weakening 
to 5%, whence the differential rotation has developed again. We discuss these simulations 
in more details in § [71 



6.2. Angular momentum balance 



The differential rotation discussed in §|1] results from angular momentum redistribution 
within the shell of intense convection realized in our simulations. In purely hydrodynamical 
models the differential rotation profile is the result of a competition among viscous torques, 
Reynolds stresses and meridional circulation. By adopting stress-free boundary conditions 
at the top and bottom boundaries of our simulations, no net torque is applied to these 
rotating spherical shells of convection. Thus total angular momentum within our simulations 
is conserved, and we can examine the ni anne r in which it is r edistr ibuted, following the 
approach developed in lEUiott et al.l (120001 ) and iBrun fc Toomrd ( 120021 ) . 



The temporal derivative of angular momentum averaged over longitude can be written 
as a divergence: 

dC _ 1 djr^J^r) ^ 1 d{sm9J^e) 
dt dr r sin 6 89 



(9) 
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with Tr and Te being respectively the radial and latitudinal components of the averaged 
angular momentum flux. They can be written and separated into their three contributions 
(due again to Viscous torque, Reynolds stresses and Meridional Circulation): 



pr sin 6 



+ve{ve + sin 9) 



Tr — J^r,V + ^r,R + J^rMC 



pr sin Q - z/r^ + v'^v'^ + Vr {v^ + sin i 



(10) 
(11) 



We can integrate the ^-fiux along concentric spheres of varying radius and the r-fiux along 
cones with various latitudinal inclinations as follows: 



n 



Igx= / J^g^xrsinedr X = V,R,MC. (12) 



lr,X 



/•TT 

/ Tr,xr^ sinOde X = V,R,MC. (13) 
Jo 

Since most of the simulations are statistically stationary, the temporal average {dC/dt)t 
vanishes. In Figure [13] we display the time-averaged integrated fluxes and Iq and their 
components for the models Yal, Ycl and Ya5T, in order to assess the effect of varying 
or/and Vt^ in the overall angular momentum balance. 

Turning first to the radial fluxes of angular momentum (Figs. [T3h . c and e), we see 
that viscous forces act to transport angular momentum radially inward in all cases. This 
inward transport is opposed mainly by the Reynolds stress flux. The role of the meridional 
circulation is less systematic. In cases Yal and Ya5T the MC flux is alternatively positive and 
negative following with good fidelity the multicellular profiles maintained in these cases and 
shown in Fig. [TTl In Ycl, the meridional circulation contribution is actually the strongest of 
all the simulations, with a characteristic amplitude two to three times larger than in most of 
the other cases. This simulation possesses a rather large and dominant clockwise cell in most 
of the domain which transports angular momentum inward at low latitude. This results in 
a meridional circulation flux that helps the viscous stresses to slow down the surface and 
speed up the base of the convection zone, thus opposing the Reynolds stresses, to yield a 
total radial angular momentum flux that is nearly zero, as noted in Fig. [12] by the solid thick 
lines. While the systems here are highly variable in time, by allowing the system to evolve for 
extended periods of time (typically thousands of days) and performing long time averages, 
we appear to be sensing the equilibrated states reasonably well. In going from the mildly 
turbulent flows of case Ycl or Ya5T to the complex ones of case Ycl, we see that the viscous 
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flux has dropped and that the Reynolds stresses and meridional circulations have changed 
accordingly to maintain equilibrium. The Reynolds stresses in cases Ycl and Ya5T are more 
evenly distributed over the whole depth of the domain, whereas in the more laminar case 
Yal is is mostly concentrated near the surface. 

Examining now the latitudinal transport of angular momentum (Figs. [T3b . d and /), 
we see that the effect of the Reynolds stresses in both cases is primarily to speed up the 
equator, since J^e,R is positive in the northern hemisphere and negative in the southern. 
It is opposed by the viscous fluxes, which act to speed up the poles and to enforce solid 
body rotation. The MC flux helps the Reynolds stresses by speeding up the equatorial 
region. This effect is opposite to what is found in si mulations of the present day Sun with 



its thinner shell of convection (IBrun &: Toomrd l2002l ) . In case Yal the MC flux is localized 



near the equator while in the two other cases these fluxes spread over a larger latitudinal 
band. In particular the circulation in model Ycl possesses a strong equatorward cell from 
low latitude (~ 20°) up to about 70°, which dominates the poleward cell near the surface. 
This larger contribution of the MC flux in this model, conflnes the Reynolds stresses toward 
lower latitudes. The manner in which each of the different components of J^e acts does not 
appear to vary appreciably in going from one case to another. There are thus no clear trends 
in following either path 7^1 (lower P^) or path #2 (higher as described in § l4.ip . However 
as the level of complexity is increased, we see a decrease in the magnitude of all components 

of J^0. 

We conclude that for these cases, involving either a low Py. or a strong rotational con- 
straint, the Reynolds stresses act latitudinally to speed up the equator, and radially to slow 
down the base of the convection zone by transporting angular momentum from the bottom 
of the domain toward the top. The latitudinal Reynolds stresses are opposed by viscous 
effects, whereas they are aided somewhat by the meridional circulations. On the contrary 
the radial Reynolds stresses have to act against both viscous and meridional flux transport. 

The overriding role of Reynolds stresses to achieve an equatorial acceleration explains, at 
least in part, the correlation between Afi and R'^ discussed in §|l]and shown Fig. [91 Here, R'^ 
can be seen as a measurement of the balance between viscous and Reynolds stresses effects. 
When i?g rises, Reynolds stresses are more efficient in transporting the angular momentum, 
since the flow is less viscous. 

This detailed balance gives us a picture of how the angular momentum is continuously 
transported, and which processes dominates acting to speed up or slow down certain regions. 
However, it is difficult to infer from the fluxes alone what the actual Vt proflle (either cylin- 
drical or conical as in the Sun) will be, besides stating that the equator will be fast or slow 
from the sense of the viscous fluxes (which are always down the gradient of VL). 
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6.3. Thermal wind 



It is well known that fluids under the influence of strong rotation can become quasi 2-D 
with most of th eir dynamics invariant with respect to dimension parallel to the rotation axis 
( jPedloskyl 119871 ) ■ Indeed if Coriolis forces dominate the others, it can be shown that: 



dz 







(14) 



which corresponds to a barotropic conflguration. It is a consequence of the so-called Taylor- 
Proudman theorem. However, rotating convection involves both radial and latitudinal heat 
transport, with the likelihood that latitudinal gradients in temperature and entropy may 
result within the convective zone. This implies that surfaces of mean press ure and density 
will riot coi ncide, thereby yielding baroclinic terms in the vorticity equations (jPedloskylll987l : 
Zahnlll992l ). Under sufficiently strong rotational constraints, a 'thermal wind balance' might 
be achieved in which departures of the angular velocity from being constant on cylinders 
(aligned with the rotation axis) are controlled by those baroclinic terms. Indeed, some 
mean-field approaches have invoked such a balance to obtain differential rotation profiles with 
bearing on t he solar convecti o n zon e fe. g.. iKitchatinov fc Rudigerlll995l : lRempelll2005al ). As 



discussed , jRn,.,fcTooU B and iMiesch et all J20n^^ore a balance 

effectively implies that 
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where z is parallel to the rotation axis. Thus latitudinal entropy gradients could serve to 
break the Taylor-Prou dman constraint which would oth erwise require the rotation to be 
constant on cylinders (jRempelll2005al : iMiesch et al.l 120061 ). However it is important to note 
that this constraint may also be broken by Reynolds and viscous stresses. Indeed the strong 
correlation found between AO and R'^ (cf. Fig. [9]) is a clear indication that Reynolds stresses 
are key players for the redistribution of angular momentum. To strengthen this point, we 
show below that these terms are as important as the baroclinic terms in establishing the 
differential rotation in the bulk of the convective envelope. 

Figure [H] assesses for case Yal and Yc5 the extent to which latitudinal entropy gra- 
dients serve to drive the temporal mean zonal flows seen as the differential rotation in 
Fig. [HI Figures [T4b and [T^ display the latitudinal entropy term seen on the right hand 
side of equation (IT^ . which in an exact thermal wind balance would be identical to dv^/dz. 
Figures [T^ anddUi show the difference between this baroclinic term and the actual dv^/dz, 
thus providing a measure of departures from such a balance. Within most of the entire 
convection zone thermal wind dominates in both cases, meaning that there is a close re- 
lationship between the zonal thermal wind and the differential rotation established in our 
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simulations. However, near the top of the domain we do see large departures from thermal 
wind balance. Further we note that for the rapidly rotating case Yc5, in the region near the 
tangent cylinder a strong shear layer (also clearly visible in Figure [TT] showing the meridional 
circulation) is present which is not in thermal wind balance. There viscous and Reynolds 
stresses play crucial roles, not unlike what i s seen in simula tions of rotating stellar radiation 
zones, the so-called Stewartson layers (e.g.. lRieutordll2006l ). 

Thus examining the nature of this thermal wind balance, combined with the assessment 
of fluxes of angular momentum, reveals that Reynolds, meridional flows and viscous stresses 
have a major role in establishing the differential rotation in the convection zone, with help 
coming from latitudinal thermal gradients. As seen in Figures H] & [5] the thermal wind is 
associated with weak latitudinal temperature variations. In Table |2] we quote the contrast 
of T at the base of the unstable domain. We find that ATbot varies from 2 K in the most 
laminar cases up to 15 K in the most turbulent ones. We also see that increasing Qq leads 
to larger temperature contrasts. However, as it can be seen in equation (fTSll . in order to 
retain a strong baroclinic term and thus a differential rotation profile that is not constant 
along cylinders, AS* (and the associated AT) must increase by exactly the same factor. 
But, as we can see in Table [21 this is hardly the case, with for instance case Yc5 having 
a contrast a factor of 2 to 3 larger than Ya5 instead of a factor of 5 that would exactly 
compensate the increase of the rotation rate. Thus the profiles become more cylindrical 
{dv^/dz 0) even though the flow remains somewhat in thermal wind balance. Thus, if 
rapidly rotating stars do not rotate along cylinders, they need an efficient latitudinal heat 
transport that establishes a strong entropy (temperature) contrast. Indeed, to get conical 
solar-like rotation profiles, it is not sufficient to have an established and balanced thermal 
wind (which should nevertheless develop after several dynamical times), but it is necessary 
to have a strong baroclinicity in the model. Such large gradients could be favored by the 
presence of a shear layer at the base of the convective domain, such as the solar tachocline 
(IMiesch et al.ll2006l ). We will indeed show in § 17.21 that enforcing a large entropy variation 
does lead to less cylindrical rotation profiles. 



7. Oscillating convection 

In this section, we focus on the oscillating behavior of convection, found in some of 
our simulations with rapid rotation, that we have mentioned several times through this 
paper. Models with low Prandtl numbers and high rotation rates, like Yc5, are typical of 
these vacillating convective solutions. We present here in more detail this phenomenon and 
discuss the physical processes governing it. 
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7.1. Prom localized to vacillating convection 

We have described in §[3] that equatorial convection tends to be localized for Yc models. 
Another phenomenon occurs in case Yc5: temporal modulations become really huge and 
the convection vacillates. The polar regions are not affected, convection there is uniformly 
developed without temporal variation. In the polar caps, the rms convective velocity w', in 
midlayer depth, remains around 24ms~^. However, in the equatorial area, convection, as 
well as being localized, oscillates between a state of fully-developed convection and a quiet 
state, characterized by very slow motions during which the convection is "asleep". In the 
equatorial band v' fluctuates between 3 to 45 ms~^. 

These oscillations are clearly not a transitional state before a statistically stationary 
state, as those discussed in § 12.41 and shown Fig. [2l The Yc5 simulation is relaxed, but 
oscillates. T his kind of behavior h as already been observed in Boussinesq simulations of the 
geodynamo ( Grote &: Busse 200 ih. and is observed here for the first time in simulations of 



deep stellar convection (see also lBrown et al.ll2007l ). 



Figure [15] shows convective patterns at two different epochs in the equatorial plane of 
the star, illustrating clearly the bimodal (high-low) state of convection. During convective 
bursts (Fig. [T5] left), convection is well developed in the equatorial plane, even though it 
is developed only on one half of the area, similar to what we see in Yc2 case. During 
quiet periods (Fig. [15] right), convection is practically fully stopped everywhere around the 
equator. We can follow these oscillations in the kinetic energy (KE) evolution (cf. Fig. [T6|l . 
This variability is almost periodic and the period is around 600 days. The averaged enthalpy 
flux (Eq. [7] , not plotted here) shows also oscillations very similar to CKE. Convective bursts 
last around 150 days, and CKE is quite fiat between each burst. The variation of differential 
rotation is different. DRKE varies like "sawteeth". It begins to increase quickly during the 
burst, so the growth lasts for ~ 150 days. Then this phase is followed by a slow decaying 
period of ~ 450 days, before a new cycle starts. As a consequence, we clearly see a phase 
difference between the maximums of DRKE and CKE. 

Left panel of Fig. [8^ allows us to see how the DRKE variation is distributed in latitude. 
We notice that, for latitudes greater than 50°, Vl is approximately uniform and does not 
change with time. Cuts of VL at latitudes 45, 60 and 75° do not show sensible differences 
between a maximum and a minimum of DRKE. This is consistent with the fact that in these 
polar regions convection is practically not modulated. Variations occur mainly around the 
equatorial band in the upper layers. It induces variations in AVt of almost 15nHz. The main 
change occurs in the slope of VL profile along the radius, mainly at latitude 0° and 15°. Such 
variations imply also variations of the radial shear in the equatorial area. 
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To investigate further, we have performed the more turbulent models Yc5T and Yd5. 
These model runs need higher resolution and thus are very expensive in terms of computing 
time. Yc5T is a model with the same Prandtl number as Yc5 but with a higher Reynolds 
number. It presents the same behavior with a slightly longer period (~ 680 days): DRKE- 
decrease phase lasts for ~ 450 days, as in case Yc5, and its bursts, slightly longer, going 
on for about 230 days typically. The persistence of the oscillating phenomenon in this more 
turbulent simulations indicates that it is not due to marginally stable conditions. All of these 
simulations are clearly supercritical. Another extra run, Yd5, as turbulent as Yc5T but with 
Pr = 1/8 shows once again oscillations. However, the period reaches ~ 1060 days, the burst 
phase length is roughly the same as for Yc5T case, but decrease duration is of ~ 800 days. Pr 
seems to influence inversely the decrease time, whereas turbulent level seems to govern the 
burst duration and so DRKE-increase length. We also notice that simulations with higher 
Pr, like Ya5 or Ya5T, do not show vacillating convection. This is because these cases have 
lower Ta (see Table [1]). The Taylor number appears to be a key parameter governing the 
oscillations; they are triggered for all simulations with a Taylor number above a threshold 
of Ta ~ 10^. Below this threshold, no models oscillate. 



Using all these results, we can build a sketch of the oscillating phenomenon. When the 
convection develops, the Reynolds stresses (v^Vj) become stronger, driving a larger differential 
rotation. In the equatorial zone, the radial shear increases (Fig. [8^). The shear is so strong 
that it destroys the coherence of convective cells (visible at midlayer depth in Fig. [T5|) . 
killing the convection. Thus (v^Vj) terms practically vanish; this implies that the differential 
rotation decreases by viscous dissipation. When the shear is sufficiently low in the equatorial 
region, convection can amplify again, and a new cycle begins. 

During quiet phases, the part of energy which is not evacuated by convection is piled 
into internal energy. During bursts, the kinetic energy is pumped from the reservoir of 
internal energy, which is on several orders of magnitude larger. With magnetic effects, we 
can imagin e very complicated s i tuations. Oscillation p roperties can drasticly change, even 



disappear (iGrote &: Bussell200ll : iMorin fc Dormyll2004l ). 



A way to verify this scenario is to follow, during the evolution, the fluxes of angular 
momentum defined in § 16. 2[ Figure [T7] shows the evolution of each integrated flux Ig along 
few oscillation cycles. We can see that near the poles — i.e., inside the tangent cylinder — 
momentum flux and all of its components are roughly null. We clearly see as expected that 
during bursts (there are three of them on the plots) the contribution of Reynolds stresses 
dominates and drags angular momentum toward the equator. Viscous effects, rather constant 
with time, act to decelerate the equator when convection has damped out, in agreement 
with the scenario described above. Effects of the meridional circulation are harder to assess 
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because it fluctuates more, especially during the bursts where the flux is rather large but its 
sign unceasingly changes. During quiet states, it helps the viscous dissipation to decelerate 
the equator, but compensates viscous effects at middle latitudes. This is a consequence of 
the multicellular topology of the meridional circulation (cf. Fig. [TT]) . 



7.2. Influence of boundary conditions: thermal forcing 

Rotation profiles of these rapidly rotating simulations are very cylindrical, like those 
predicted by the first solar simulations before they were negated by helioseismic observa- 
tions. Up to now, there has been no observational inference of the inner rotation profile 
for other stars. Thus, such cylindrical profiles could reflect reality, but we cannot exclude 
the possibility that our simulations suffer problems similar to those of the first solar simula- 
tions. In this section, we want to see how the dynamics is modified when the "cylindricity" 
of rotation profiles is reduced. Indeed, if the rotation profiles of thermally forced models 
are more conical than those of previous unforced models, the shear in the equatorial plane 
would have to be reduced, along with the oscillation process. Profiles of the Yc5 model are 
almost cylindrical, due to the high rotation rate. Coriolis forces dominate the others and the 
equation f|T^ is almost verified. The Taylor-Proudman theorem can be broken if a baroclinic 
term exists, as defined by the equation f[T^ . However, as seen in § 16.31 this term has reduced 
in strength in going from Ycl to Yc5. Thus a potential thermal wind can lead to more 
conical profiles (see discussion § 16. 3p . The origin of such a therrnal wind could be a strongl y 



sheared region, like the tachocline observed in the Sun (lRempelll2005al : iMiesch et al.ll2006l ). 



We have imposed such an entropy gradient at the base of the convective shell in a twin 
of Yc5 simulation: Yc5S. The latitudinal entropy profile imposed at the bottom of the shell 
is in the form 

02^2° + ^4^°, (16) 



'S'bot 



Cp 

with 02 = 1.14 X 10~^ and = — 1.8 x 10^^. In the present configuration, compared to 



solar case (IMiesch et al.ll2006l ). profiles are strongly cylindrical and the rotation speed is five 
times higher. So we need to impose a strong gradient to expect a sensible effect. In this 
simulation the contrast of temperature between poles and equator is around 60 K, that is 5 
times larger than the Sun, 9 times larger than Ycl case, and 3.5 larger than the initial Yc5 
model (see also Table [2]). 

Figure [TS] shows the temporal evolution of kinetic energy. After a cycle similar to 
those shown by Yc5, the entropy has diffused and we reached a phase characterized by a 
reduction of the oscillation amphtudes, as we had expected. Figure [8]/ shows this profile at 
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t ~ 2000 days. The rotation is less cylindrical, even less than Ycl case. As a consequence, 
the contrast of differential rotation along radius decreases in the equatorial zone, which is 
clearly seen by comparing the slopes of curves at 0° in the right panels of Fig. [H^ and[H]/. 
As a function of latitude, the contrast of rotation increases and the profile becomes purely 
monotonic: the decrease of rotation continues at high latitudes — i.e., inside the tangent 
cylinder — and does not stop around 50° as in Yc5. However, after t ~ 2600 days, our model 
achieves a rather chaotic new state (Fig. [T8|) . There are strong variations in the convective 
state again, but no more simple periodicity. Indeed, differential rotation has continued to 
increase in latitude, but in radius also, making possible a sufficiently strong shear again. 

As expected, imposing an entropy gradient at the base of the convective zone could lead 
to a more conical rotation profile, modifying the dynamics and reducing the oscillations. The 
thermal c ontras t need ed can be caused by a strong shear due to a thinner tachocline as was 



shown by iBrunI (120071 . eq. [11]). Indeed, in order to get a more conical profile for the same 
mean radial jump in angular velocity in the tachocline, one needs at fixed rotation rate to have 
a larger latitudinal entropy variation or equivalently a thinner tachocline. Indications of the 
nature of profiles should be obtained from asteroseismology, in the next years. Even although 
the quality of asteroseismic data will not be comparable to that available fo r helioseismology. 



2003; 


Ballot. Garcia. & Lambert 


2006), and the 


inferred for rapid rotators i 


Gizon & Solanki 


2004) 



8. Conclusions and perspectives 

Our simplified 3-D numerical simulations of the turbulent convective envelope in young 
solar-like stars have served to catch some aspects of the very complex dynamics which occurs 
in such thick shells and show us how convection and rotation are intricate, even if we are, for 
sure, far from a true young star. We have especially seen how these simulations are sensitive 
to the Taylor number T^. 

The main important result of this work concerns the behavior of mean fiows according to 
the rotation speed of the star. We find that the contrast of the differential rotation decreases 
only slightly, whereas the intensity of the meridional circulation clearly decays. High 
make our rotation profiles of rapidly rotating envelopes highly cylindrical. It is possible to 
make them more conical — as exhibited by the true Sun — by imposing at the bottom of 
the convective zone a thermal forcing in latitude, such as may arise from a sharp stellar 
tachocline. 
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While global convective patterns are quite similar to those seen in solar simulations, 
some clear differences have appeared. For lower hot spots emerg e (see Ycl). Such s pots 
are also seen in a thin solar envelope, but with higher rotation rates (IBrown et al.ll2007l ). At 
higher rotation rates, hot spots change into localized convection, where convective motions 
are preferentially developed in a particular range of longitudes in the equatorial zone. Oscil- 
lations in convection appear in models characterized by high Tq. The shell thickness has to 
play a role, by favoring the appearance of a strong-shear layer in the n iiddle of the shell. Th e 
situation is similar to those also observed in geophysics simulations (iGrote fc Bussell200ll ). 
The mechanism of oscillations is well understood, especially the role of the Reynolds stresses 
that drive a strong differential rotation which sparks off the shear. If this kind of vacillating 
and/or localized deep convection occurs in real stars, their surfaces should be affected. It 
could be seen in observations — of luminosity, line widths... — as temporal variations, aris- 
ing from changing in the convective state and/or from the rotation of different longitudes 
into the line of sight. Young stars, especially T Tauri stars, are known to be variable objects, 
but the luminosity fluctuations are linked to magnetic field (chromospheric or coronal activ- 
ity, photospheric spots...) or, in some cases, accretion phenomena. In the present situation, 
it would be specious to make any comparison of the deep-convection behavior seen in our 
most rapidly rotating simulations, with observed stellar luminosity curves (dominated by 
surface effects). Nevertheless, if we focus on the differential rotation variability seen in sim- 
ulations, we note that significant secular variations of differential rotatio n contrast have also 
been observed in young K dwarfs with high rotation rate like AB Dor (jCameron &: Donati 
2OO2I : iDonati. Cameron. &: Petitll2003al ). However, these observed variations can also have a 
magnetic origin. 

The next step of this work would be to see the effect of these specific properties both 
of convection and of mean flows, on the dynamo processes which can occur in these stars. 
These purely hydrodynamical simulations have exhibited obvious differences with the solar 
ones, that let forshadow also differences with the Sun, when the magnetic field is included. 
Observations tend to show that magnetic properties of young stars are peculiar. All of them 
show a magnetic activity which can be tracked by their X-ray emissions. Several issues have 
been addressed about the signification of correlations — or the lack of correlation — between 



rotat i on rate a.nd ma gnetic activity. Is it just an extension of the saturation (IStauffer et al. 



1994 iRandichI 119971 ) and "supersaturation" phenomena (jProsser et al.l Il996 



James et al. 



2OOOI ) observed for main sequence stars? Or is it the signature o f a turbulent dynamo 



distributed throughout t he deep convectiv e zone (Durnev et al.l 19931 ) instead of a classic a- 
f2 dynamo (Parker 1993h. as suggested by Feigelson et al. ( 2003 )? Zeeman-Doppler imaging 
of young fast rotators (IDonati et al.ll2003bl ) shows large-scale azimuthal magnetic structures, 
that seems to plead for such dynamo mechanisms. Such fascinating issues encourage us to 
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pursue simulations of young stars by including magnetic field, and searching for clues to the 
answers of these, and other, questions. 

First we acknowledge useful discussions with J. Toomre and B. Brown. We also wish 
to thank the referee, P. Oilman, for helping to clarify the paper. The simulations were 
performed using the computer facilities of CEA/CCRT. J.B. is grateful to the developers of 
the ASH code for letting him use it for this study. 
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Table 1. Parameters and characteristics of models 



Model 


n/nQ 






u 








K 




Pr 








Ta 


Roc 


Ro 


Re/R'e 




Yal . . . 


1 


3, 


,7 


X 


1012 


9, 


.3 


X 


10" 


4 


1, 


.7 


X 


iqs 


1.8 X 10'' 


0.15 


0.013 


22/18 


65, 128, 256 


Ya2 . . . 


2 


3, 


,7 


X 


1012 


9, 


.3 


X 


10" 


4 


3, 


,8 


X 


iqs 


7.2 X 106 


0.11 


0.0053 


18/15 


65,128,256 


Ya5 . . . 


5 


3. 


,7 


X 


1012 


9, 


.3 


X 


10" 


4 


1, 


,1 


X 


id" 


4.5 X 10'' 


0.08 


0.0016 


13/11 


65,128,256 


Ya5T . 


5 


9, 


,1 


X 


10" 


2, 


,3 


X 


10" 


4 


1, 


,2 


X 


10^ 


7.5 X 10** 


0.06 


0.0020 


119/54 


129, 256,512 


Ybl . . 


1 


1, 


,9 


X 


1012 


1, 


.9 


X 


1012 


1 


4, 


,0 


X 


IQS 


7.2 X 10" 


0.23 


0.014 


83/38 


65, 128, 256 


Yb2 . . 


2 


1, 


,9 


X 


1012 


1, 


,9 


X 


1012 


1 


8, 


,6 


X 


IQS 


2.9 X 10^ 


0.17 


0.0059 


72/31 


65,128,256 


Yb5 . . 


5 


1, 


,9 


X 


1012 


1, 


,9 


X 


1012 


1 


2, 


,4 


X 


ID'S 


1.8 X 10** 


0.11 


0.0020 


46/27 


65,128,256 


Yb5T. 


5 


9, 


,3 


X 


10" 


9, 


,3 


X 


10" 


1 


7, 


.4 


X 


10*^ 


7.5 X 10* 


0.10 


0.0020 


155/55 


129, 256,512 


Ycl . . . 


1 


4, 


,6 


X 


10" 


1, 


.8 


X 


1012 


1/4 


2, 


,7 


X 


10*5 


1.2 X 10* 


0.29 


0.020 


628/224 


129, 256,512 


Yc2. . . 


2 


4. 


,6 


X 


10" 


1, 


,8 


X 


1012 


1/4 


5, 


.0 


X 


10^ 


4.8 X 10* 


0.20 


0.0084 


557/180 


129, 256,512 


Yc5. . . 


5 


4, 


,6 


X 


10" 


1, 


,8 


X 


1012 


1/4 


1, 


.1 


X 


10^ 


3.0 X 10'' 


0.12 


0.0024=* 


490/128=* 


129, 256,512 


Yc5St 


5 


4, 


,6 


X 


10" 


1, 


.8 


X 


1012 


1/4 


1, 


,1 


X 


10^ 


2.9 X 10^ 


0.12 


0.0025^ 


728/138=* 


129, 256,512 


Yc5T. 


5 


2, 


,3 


X 


10" 


9, 


.3 


X 


10" 


1/4 


3, 


,3 


X 


10^ 


1.2 X 10"' 


0.11 


0.0024^ 


1253/256=* 


161,512, 1024 


Yd5 . . 


5 


2. 


,3 


X 


10" 


1, 


,8 


X 


1012 


1/8 


2, 


.1 


X 


10^ 


1.2 X 10"' 


0.12 


0.0032=* 


1297/345=* 


161,512, 1024 



Note. — All simulations have an inner radius rj, = 4.2 X lOl" cm and an outer radius rt = 7.3 X lOl" cm, what is to say a 
thickness of the computational domain L = 3.1 X IQi" cm. The eddy viscosity u and conductivity k at mid-depth are quoted in 
cm2 s~l. Here evaluated at mid-layer depth are the temporal average of the Prandtl number Pr = i'/k, the Rayleigh number 
Ra = -{dS /dr){dp/dS)gL* /(pvK), the Taylor number Ta = 402^4/^^2^ ^j^g convective Rossby number Roo = [-Ra/(TaP^)]l/2, 
the rms Rossby number Ro = v' /{2QoL), and the rms Reynolds numbers Re = vL/u and R'^ = v'L/u, where v and v' are the 
rms velocity and rms convective velocity at middepth (cf. Table [2J. The number of radial, latitudinal and longitudinal mesh 
points are Nr, Ng, N^. 

t Model with a thermal forcing (see 17.21 1. 

=*Due to vacillating behavior of convection in these models, rms velocities fluctuate noticeably, thus Ro and Re vary a lot 
according to the convective level. 



Table 2. Representative velocities, energy contents, differential rotation, and temperature contrast 
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Note. — Temporal averages of the rms velocity v, of the rms components Cr, fe, f0, and of fluctuating velocities v' and v'^ (axisymmetric 
components are removed) are estimated at midlaycr depth, all expressed in units of ms~^. ar = v'^'^/v''^ is the anisotropy index. Also listed are 
the time average over the complete volume of the total kinetic energy KE and that associated with the (axisymmetric) differential rotation DRKE, 
the (axisymmetric) meridional circulation MCKE, and the non-axisymmetric convection CKE, all in units of ergcm"^. The latitudinal contrasts 
between 0° and 60° of angular frequencies AQ (resp. temperatures T), quoted in nHz (resp. K), axe computed at the top (resp. bottom) of the 
domain. 
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Fig. 1. — Thermodynamic mean radial stratification in the outer part of a young 1-Mq 
star, within both the 1-D stellar model Y and 3-D hydrodynamic case Ybl. The mean 
temperature T, density p and pressure P are plotted as curves for the model Y and with 
symbols at their mesh locations for our case Ybl. The base of the outer convective zone 
(bcz) is indicated by the dotted line. Variables are normalized to the value they have in 
model Ybl at the bottom of the computational domain: Ti, — 3.4 x 10^ K, pi, — 1.9 gcm~^, 
and Pb — 9.1 x 10^^ dynecm"^. 
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Fig. 2. — Evolution of the enthalpy flux and of the kinetic energy density (KE), averaged 
over the full domain from the quiescent state to a statistically stationary state {solid lines). 
The evolution of the three components of KE are also plotted: the kinetic energy in convective 
motions {dashed line), in differential rotation {dotted line), and in meridional circulation 
{dash-dotted line). 
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Fig. 3. — Radial transport of energy in case Ycl achieved by the fluxes Le„, L^d, Lfee, l^ed-, 
L^, and their total L^t, all normalized by the stellar luminosity L*. 
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Fig. 4. — Radial velocity Vr and temperature fluctuation T flelds near the surface (0.93 i?*) 
and at mid-depth (0.75 at a given instant for the three models Yal, Ya5, and YaST. Maps 
are shown in MoUweide projection. Dashed horizontal line designates equator. Meridians 
are plotted with dotted lines every 45° and parallels every 30°. Dotted ellipse shows the 
relative position of the stellar surface. Minimum and maximum amplitudes for each fleld are 
indicated next to each panel. 
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Fig. 5. — Fields Vr and T near the surface and at mid-depth, for the three models Ycl, Yc2, 
and Yc5. See caption of Fig. HI 
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Fig. 6. — Entropy gradient dS/dr for models Yal {solid line), Ya2 {dashed line), Ya5 {dash- 
dotted line), and Ya5T {dotted line). 
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Fig. 7. — (0, t)-diagrams showing the temporal evolution of the Vr profile along the equator 
{9c = 90°), close to the top of the shell (rc = 0.93 R^^) for models Yal (a, 6) and Ycl 
{c,d). On the left (a, c), diagrams are plotted in the corotating frame. On the right {b,d), 
the diagrams are plotted in the "shifted" {4>s,t) coordinates, to remove the local effect of 
differential rotation (see text). 
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Fig. 8. — Temporal and azimuthal averages of angular velocities achieved for six simulations: 
Yal, Ybl, Ycl, Yc2, Yc5, and Yc5S. For each panel a-f, a contour plot of Q/27r in the 
meridional plane is shown on the left, and on the right, radial cuts of the same quantity are 
plotted at six fixed latitudes (specified on panel c) , averaged on both hemispheres. The color 
scale of each map is independent and is indicated nearby the plot. Thick lines on contour 
plots correspond to i7 = Qo- The scales of the right plots are the same for panels a-d, and 
also for panels e and /. Since Yc5 is an oscillating model, the e right panel shows radial 
cuts of fl/2TT at two different moments: when the differential rotation reaches 1) a maximum 
[solid lines) and 2) a minimum [dashed lines). 
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Fig. 9. — Differential rotation contrasts AQ according to convective Reynolds numbers R'^. 
Plus, diamonds, crosses, and square correspond to Ya, Yb, Yc, and Yd models respectively. 
T indexes indicate models Ya5T, Yb5T, and Yc5T. 
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Fig. 10. — Differential rotation contrasts AQ according to rotation rates Qo- The sense of 
symbols is the same as in Fig. [9l Dotted lines correspond to the scaling law deduced by 
logarithmic regression for each series. Values of fitted slopes are printed out on the right. 
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Fig. 11. — Meridional circulation in Yal, Ycl, and Yc5. Maps show the mean mass flux cir- 
culation in the azimuthal plane. White lines over black background correspond to clockwise 
circulations, and black lines over white background to counterclockwise ones. The color scale 
is not linear to make visible even weak flows. Mean latitudinal velocity vg profiles at the top 
of the shell are also plotted. Maps and profiles are obtained by averaging over longitude and 
time. 
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Fig. 12. — Surface velocity w^max characterizing the meridional circulation intensity, plotted 
as function of the rotation rate flo- See caption of Fig. [TDl 
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Fig. 13. — Time average of the latitudinal integral of the radial angular momentum flux 
Ir (left), and of the radial integral of the latitudinal angular momentum flux Iq (right), for 
models Yal (top), Ycl (middle) and Ya5T {bottom). The fluxes have been decomposed into 
their viscous (dotted lines), Reynolds stress (dot-dashed lines), and meridional circulation 
components (dot-dashed lines). The total fluxes are plotted with solid lines. 
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Fig. 14. — Temporal and azimuthal averages for models Yal (a, b) and Yc5 (c, d) showing the 
baroclinic term in the meridional force balance — defined eq. f|T5|) — (a, c), and the difference 
of this term with dVfp/dz {b, d). 




Fig. 15. — Radial velocity field in the equatorial plane for model Yc5 during a maximum 
{left) and a minimum [right] of convective activity. 
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Fig. 16. — Temporal evolution of kinetic energy content for the Yc5 case. KE {solid line) 
is decomposed into its components: DRKE {dotted line), CKE {dashed line), and MCKE 
{dot-dashed line). 
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Fig. 17. — Evolution as a function of time of 1$, the angular-momentum flux in 6'-direction, 
averaged on radius. Ig {top left) is decomposed into the contributions of Reynolds stresses 
h,R {top right), meridional circulation Ig^MC {bottom left), and viscous torque Igy {bottom 
right). White/red (black/green) indicate southward (northward) transfer. 
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Fig. 18. — Temporal evolution of kinetic energy content for the Yc5S case. See caption of 
Fig.[ISl 



